function manifold_analysis_visualization_functions(datain, fflag)
    %%% global parameters %%%
    Fsi = 20;
    clrs = [0, 0.4470, 0.7410; 0.8500, 0.3250, 0.0980; 0.9290, 0.6940, 0.1250; 0.4940, 0.1840, 0.5560; 0.4660, 0.6740, 0.1880; 0.3010, 0.7450, 0.9330; 0.6350, 0.0780, 0.1840; 0.25, 0.25, 0.25];

    switch fflag
        case 'sync'
            wid = 0.5;
            clf
            hold on
            for i = 1: length(a)
                tmp = a{i}(2, :);
                plot(i * ones(size(tmp)) + wid * rand(size(tmp)) - wid / 2, tmp, '.k')
                errorbar(i, mean(tmp), std(tmp) / sqrt(length(tmp)) / 2, '*r')
            end
            ylim([0, 0.01])
            
        case 'cluster_major_property'
            dt = datain{1};
            whwpp = dt.whwpp;
            whwpn = dt.whwpn;
            whnwpp = dt.whnwpp;
            whnwpn = dt.whnwpn;
            nwhwpp = dt.nwhwpp;
            nwhwpn = dt.nwhwpn;
            nwhnwpp = dt.nwhnwpp;
            nwhnwpn = dt.nwhnwpn;
            [nc, nw] = size(whwpp);
            
            figure(gcf)
            clf
            count = 1;
            for i = 1: nc
                for j = 1: nw
                    subplot(nc, nw, count)
                    hold on
                    xrg = [];
                    mn = [];
                    stdd = [];
                    for ii = 1: length(whwpp{i, j})
                        if ~isempty(whwpp{i, j}{ii})
%                             plot(ii, whwpp{i, j}{ii}(:, 2), '.', 'color', clrs(1, :))
%                             plot(ii, whnwpp{i, j}{ii}(:, 2), '.', 'color', clrs(2, :))
%                             plot(ii, nwhwpp{i, j}{ii}(:, 2), '.', 'color', clrs(3, :))
%                             plot(ii, nwhnwpp{i, j}{ii}(:, 2), '.', 'color', clrs(4, :))
                            xrg = [xrg, ii];
                            mn = [mn, [mean(whwpp{i, j}{ii}(:, 2)); mean(whnwpp{i, j}{ii}(:, 2)); mean(nwhwpp{i, j}{ii}(:, 2)); mean(nwhnwpp{i, j}{ii}(:, 2))]];
                            stdd = [stdd, [std(whwpp{i, j}{ii}(:, 2)) / (2 * sqrt(size(whwpp{i, j}{ii}, 1))); std(whnwpp{i, j}{ii}(:, 2)) / (2 * sqrt(size(whnwpp{i, j}{ii}, 1))); std(nwhwpp{i, j}{ii}(:, 2)) / (2 * sqrt(size(nwhwpp{i, j}{ii}, 1))); std(nwhnwpp{i, j}{ii}(:, 2)) / (2 * sqrt(size(nwhnwpp{i, j}{ii}, 1)))]];
                        end
                    end
                    h1 = errorbar(xrg, mn(1, :), stdd(1, :), '*-', 'color', clrs(1, :));
                    h2 = errorbar(xrg, mn(2, :), stdd(2, :), '*-', 'color', clrs(2, :));
                    h3 = errorbar(xrg, mn(3, :), stdd(3, :), '*-', 'color', clrs(3, :));
                    h4 = errorbar(xrg, mn(4, :), stdd(4, :), '*-', 'color', clrs(4, :));
                    legend([h1, h2, h3, h4], {'whwp', 'whnwp', 'nwhwp', 'nwhnwp'}, 'location', 'best')
                    count = count + 1;
                end
            end
            
        case 'cluster_major_property_2d_all_time'
            dt = datain{1};
            whwpp = dt.whwpp;
            whwpn = dt.whwpn;
            whnwpp = dt.whnwpp;
            whnwpn = dt.whnwpn;
            nwhwpp = dt.nwhwpp;
            nwhwpn = dt.nwhwpn;
            nwhnwpp = dt.nwhnwpp;
            nwhnwpn = dt.nwhnwpn;
            [nc, nw] = size(whwpp);
            
            figure(gcf)
            clf
            count = 1;
            for i = 1: nc
                for j = 1: nw
                    subplot(nc, nw, count)
                    hold on
                    xrg = [];
                    mn1 = [];
                    mn2 = [];
                    for ii = 1: length(whwpp{i, j})
                        if ~isempty(whwpp{i, j}{ii})
%                             plot(ii, whwpp{i, j}{ii}(:, 2), '.', 'color', clrs(1, :))
%                             plot(ii, whnwpp{i, j}{ii}(:, 2), '.', 'color', clrs(2, :))
%                             plot(ii, nwhwpp{i, j}{ii}(:, 2), '.', 'color', clrs(3, :))
%                             plot(ii, nwhnwpp{i, j}{ii}(:, 2), '.', 'color', clrs(4, :))
                            xrg = [xrg, ii];
                            mn1 = [mn1, [mean(whwpp{i, j}{ii}(:, 1)); mean(whnwpp{i, j}{ii}(:, 1)); mean(nwhwpp{i, j}{ii}(:, 1)); mean(nwhnwpp{i, j}{ii}(:, 1))]];
                            mn2 = [mn2, [mean(whwpp{i, j}{ii}(:, 2)); mean(whnwpp{i, j}{ii}(:, 2)); mean(nwhwpp{i, j}{ii}(:, 2)); mean(nwhnwpp{i, j}{ii}(:, 2))]];
                        end
                    end
                    h1 = plot(mn1(1, :), mn2(1, :), 'color', clrs(1, :));
                    h2 = plot(mn1(2, :), mn2(2, :), 'color', clrs(2, :));
                    h3 = plot(mn1(3, :), mn2(3, :), 'color', clrs(3, :));
                    h4 = plot(mn1(4, :), mn2(4, :), 'color', clrs(4, :));
                    legend([h1, h2, h3, h4], {'whwp', 'whnwp', 'nwhwp', 'nwhnwp'}, 'location', 'best')
                    count = count + 1;
                end
            end

        case 'cluster_major_property_2d_after_errorbar'
            dt = datain{1};
            whwpp = dt.whwpp;
            whwpn = dt.whwpn;
            whnwpp = dt.whnwpp;
            whnwpn = dt.whnwpn;
            nwhwpp = dt.nwhwpp;
            nwhwpn = dt.nwhwpn;
            nwhnwpp = dt.nwhnwpp;
            nwhnwpn = dt.nwhnwpn;
            [nc, nw] = size(whwpp);
            tstep = [4, 1];
            
            figure(gcf)
            clf
            count = 1;
            for i = 1: nc
                for j = 1: nw
                    subplot(nc, nw, count)
                    hold on
                    xrg = [];
                    x1 = [];
                    x2 = [];
                    x3 = [];
                    x4 = [];
                    y1 = [];
                    y2 = [];
                    y3 = [];
                    y4 = [];
                    for ii = 2: tstep(i) + 1
                        if ~isempty(whwpp{i, j}{ii})
%                             plot(ii, whwpp{i, j}{ii}(:, 2), '.', 'color', clrs(1, :))
%                             plot(ii, whnwpp{i, j}{ii}(:, 2), '.', 'color', clrs(2, :))
%                             plot(ii, nwhwpp{i, j}{ii}(:, 2), '.', 'color', clrs(3, :))
%                             plot(ii, nwhnwpp{i, j}{ii}(:, 2), '.', 'color', clrs(4, :))
                            xrg = [xrg, ii];
                            x1 = [x1; whwpp{i, j}{ii}(:, 1)];
                            x2 = [x2; whnwpp{i, j}{ii}(:, 1)];
                            x3 = [x3; nwhwpp{i, j}{ii}(:, 1)];
                            x4 = [x4; nwhnwpp{i, j}{ii}(:, 1)];
                            y1 = [y1; whwpp{i, j}{ii}(:, 2)];
                            y2 = [y2; whnwpp{i, j}{ii}(:, 2)];
                            y3 = [y3; nwhwpp{i, j}{ii}(:, 2)];
                            y4 = [y4; nwhnwpp{i, j}{ii}(:, 2)];
                        end
                    end
                    xmn1 = mean(x1);
                    xmn2 = mean(x2);
                    xmn3 = mean(x3);
                    xmn4 = mean(x4);
                    ymn1 = mean(y1);
                    ymn2 = mean(y2);
                    ymn3 = mean(y3);
                    ymn4 = mean(y4);
                    xst1 = std(x1) / (2 * sqrt(length(x1)));
                    xst2 = std(x2) / (2 * sqrt(length(x2)));
                    xst3 = std(x3) / (2 * sqrt(length(x3)));
                    xst4 = std(x4) / (2 * sqrt(length(x4)));
                    yst1 = std(y1) / (2 * sqrt(length(y1)));
                    yst2 = std(y2) / (2 * sqrt(length(y1)));
                    yst3 = std(y3) / (2 * sqrt(length(y1)));
                    yst4 = std(y4) / (2 * sqrt(length(y1)));
                    h1 = errorbar(xmn1, ymn1, yst1, yst1, xst1, xst1, 'color', clrs(1, :));
                    h2 = errorbar(xmn2, ymn2, yst2, yst2, xst2, xst2, 'color', clrs(2, :));
                    h3 = errorbar(xmn3, ymn3, yst3, yst3, xst3, xst3, 'color', clrs(3, :));
                    h4 = errorbar(xmn4, ymn4, yst4, yst4, xst4, xst4, 'color', clrs(4, :));
                    legend([h1, h2, h3, h4], {'whwp', 'whnwp', 'nwhwp', 'nwhnwp'}, 'location', 'best')
                    count = count + 1;
                end
            end

        case 'cluster_major_property_2d_combined'
            dt = datain{1};
            whwpp = dt.whwpp;
            whwpn = dt.whwpn;
            whnwpp = dt.whnwpp;
            whnwpn = dt.whnwpn;
            nwhwpp = dt.nwhwpp;
            nwhwpn = dt.nwhwpn;
            nwhnwpp = dt.nwhnwpp;
            nwhnwpn = dt.nwhnwpn;
            [nc, nw] = size(whwpp);
            tstep = [4, 1];
            mks = {'o', 's', '*', 'h'};
            lgds = {{'heat whwp', 'heat whnwp', 'heat nwhwp', 'heat nwhnwp', 'heat no whwp', 'heat no whnwp', 'heat no nwhwp', 'heat no nwhnwp'}, {'vf whwp', 'vf whnwp', 'vf nwhwp', 'vf nwhnwp', 'vf no whwp', 'vf no whnwp', 'vf no nwhwp', 'vf no nwhnwp'}};
            ttl = {'heat', 'vf'};
            
            figure(gcf)
            clf
            for i = 1: nc
                count = 1;
                hs = [];
                subplot(2, 1, i)
                hold on
                x1a = [];
                x2a = [];
                x3a = [];
                x4a = [];
                y1a = [];
                y2a = [];
                y3a = [];
                y4a = [];
                for j = 1: nw
                    xrg = [];
                    x1 = [];
                    x2 = [];
                    x3 = [];
                    x4 = [];
                    y1 = [];
                    y2 = [];
                    y3 = [];
                    y4 = [];
                    for ii = 2: tstep(i) + 1
                        if ~isempty(whwpp{i, j}{ii})
%                             plot(ii, whwpp{i, j}{ii}(:, 2), '.', 'color', clrs(1, :))
%                             plot(ii, whnwpp{i, j}{ii}(:, 2), '.', 'color', clrs(2, :))
%                             plot(ii, nwhwpp{i, j}{ii}(:, 2), '.', 'color', clrs(3, :))
%                             plot(ii, nwhnwpp{i, j}{ii}(:, 2), '.', 'color', clrs(4, :))
                            xrg = [xrg, ii];
                            x1 = [x1; whwpp{i, j}{ii}(:, 1)];
                            x2 = [x2; whnwpp{i, j}{ii}(:, 1)];
                            x3 = [x3; nwhwpp{i, j}{ii}(:, 1)];
                            x4 = [x4; nwhnwpp{i, j}{ii}(:, 1)];
                            y1 = [y1; whwpp{i, j}{ii}(:, 2)];
                            y2 = [y2; whnwpp{i, j}{ii}(:, 2)];
                            y3 = [y3; nwhwpp{i, j}{ii}(:, 2)];
                            y4 = [y4; nwhnwpp{i, j}{ii}(:, 2)];
                        end
                    end
                    xmn1 = mean(x1);
                    xmn2 = mean(x2);
                    xmn3 = mean(x3);
                    xmn4 = mean(x4);
                    ymn1 = mean(y1);
                    ymn2 = mean(y2);
                    ymn3 = mean(y3);
                    ymn4 = mean(y4);
                    xst1 = std(x1) / (2 * sqrt(length(x1)));
                    xst2 = std(x2) / (2 * sqrt(length(x2)));
                    xst3 = std(x3) / (2 * sqrt(length(x3)));
                    xst4 = std(x4) / (2 * sqrt(length(x4)));
                    yst1 = std(y1) / (2 * sqrt(length(y1)));
                    yst2 = std(y2) / (2 * sqrt(length(y1)));
                    yst3 = std(y3) / (2 * sqrt(length(y1)));
                    yst4 = std(y4) / (2 * sqrt(length(y1)));
                    
                    x1a = [x1a, xmn1];
                    x2a = [x2a, xmn2];
                    x3a = [x3a, xmn3];
                    x4a = [x4a, xmn4];
                    y1a = [y1a, ymn1];
                    y2a = [y2a, ymn2];
                    y3a = [y3a, ymn3];
                    y4a = [y4a, ymn4];
                    
                    hs = [hs, errorbar(xmn1, ymn1, yst1, yst1, xst1, xst1, mks{count}, 'color', clrs(1, :))];
                    hs = [hs, errorbar(xmn2, ymn2, yst2, yst2, xst2, xst2, mks{count}, 'color', clrs(2, :))];
                    hs = [hs, errorbar(xmn3, ymn3, yst3, yst3, xst3, xst3, mks{count}, 'color', clrs(3, :))];
                    hs = [hs, errorbar(xmn4, ymn4, yst4, yst4, xst4, xst4, mks{count}, 'color', clrs(4, :))];
                    count = count + 1;
                end
                plot(x1a, y1a, 'color', clrs(1, :))
                plot(x2a, y2a, 'color', clrs(2, :))
                plot(x3a, y3a, 'color', clrs(3, :))
                plot(x4a, y4a, 'color', clrs(4, :))
                title(ttl{i})
                xlabel('Population active level')
                ylabel('Synchrony with Gaussian kernel')
                legend(hs, lgds{i}, 'location', 'best')
            end
            
        case 'trajectory_length_full_plot'
            tps = datain{1};
            [nc, nw, nm] = size(tps);
            aa = squeeze(tps(:, 1, :));
            bb = cell(nc, 4);
            for i = 1: nc
                for j = 1: 4
                    for k = 1: nm
                        tmp = aa{i, k}{j};
                        bb{i, j} = [bb{i, j}; tmp];
                    end
                end
            end
            
            for i = 1: nc
                tmp = bb(i, :);
                subplot(nc, 1, i)
                mn = cellfun(@(x) nanmean(x, 1), tmp);
                stdd = cellfun(@(x) nanstd(x, [], 1) / (2 * sqrt(size(x, 1))), tmp);
                errorbar(mn, stdd)
                ylim([0.96, 1.12])
            end
            
        case 'visualize_geodesic'
            dt = datain{1};
            [nc, nw, nm] = size(dt);
            tstep = [1, 6; 1, 2] * Fsi;
            
            figure(gcf)
            clf
            count = 1;
            for i = 2: nc
                for j = 1: nw
                    tp = cell(1, 4);
                    for k = 1: nm
                        for ii = 1: 4
                            tp{ii} = [tp{ii}; dt{i, j, k}{ii}];
                        end
                    end
                    subplot(nc - 1, nw, count)
                    hold on
                    mx = 0;
                    for ii = 1: 4
                        tmp = tp{ii};
%                         tmp = tp{ii}(:, tstep(i - 1, 1): tstep(i - 1, 2));
                        tmp = cumsum(tmp, 2);
%                         tmp = cumsum(tmp, 2) - tmp(:, 1);
%                         plot(tmp', 'color', clrs(ii, :))
                        plot(mean(tmp, 1), 'color', clrs(ii, :), 'linewidth', 3)
                        mx = max(max(mean(tmp, 1)), mx);
                    end
%                     ylim([0.8 * mx, 1.1 * mx])
                    count = count + 1;
                end
            end

        case 'trajectory_length'
            dt = datain{1};
            [nc, nw, nm] = size(dt);
            tstep = [5, 2] * Fsi;
            ttls = {'heat full', 'heat no'; 'vf full', 'vf no'};
            
            tp = cell(nc - 1, nw);
            count = 1;
            figure(gcf)
            clf
            for i = 2: nc
                for j = 1: nw
                    tp{i, j} = NaN(nm, 4);
                    for k = 1: nm
                        dtmp = dt{i, j, k};
                        for ii = 1: length(dtmp)
                            tmp = dtmp{ii};
                            if ~isempty(tmp)
                                rt = sum(tmp(:, Fsi + 1: tstep(i - 1)), 2);
                                tp{i, j}(k, ii) = mean(rt);
                            end
                        end
                        tp{i, j} = tp{i, j} ./ tp{i, j}(:, 1);
                    end
                    
                    tmp = tp{i, j};
                    subplot(nc - 1, nw, count)
                    mn = nanmean(tmp, 1);
                    stdd = nanstd(tmp, [], 1) / (2 * sqrt(size(tmp, 1)));
                    errorbar(mn, stdd)
                    title(ttls{i - 1, j})
                    count = count + 1;
                end
            end
            
        case 'trajectory_length_full'
            dt = datain{1};
            [nc, nw, nm] = size(dt);
            tstep = [5, 2] * Fsi;
            ttls = {'heat full', 'heat no'; 'vf full', 'vf no'};
            
            tp = cell(nc - 1, nw);
            count = 1;
            figure(gcf)
            clf
            for i = 2: nc
                for j = 1: nw
                    tp{i, j} = cell(1, 4);
%                     for k = 1: nm
%                         dtmp = dt{i, j, k};
%                         for ii = 1: length(dtmp)
%                             tmp = dtmp{ii};
%                             if ~isempty(tmp)
%                                 rt = sum(tmp(:, Fsi + 1: tstep(i - 1)), 2);
%                                 tp{i, j}{ii} = [tp{i, j}{ii}; rt];
%                             end
%                         end
%                         tp{i, j} = cellfun(@(x) x / mean(tp{i, j}{1}), tp{i, j}, 'uniformoutput', false);
%                     end
                    for k = 1: nm
                        dtmp = dt{i, j, k};
                        rtt = cell(1, length(dtmp));
                        for ii = 1: length(dtmp)
                            tmp = dtmp{ii};
                            if ~isempty(tmp)
                                rt = sum(tmp(:, Fsi + 1: tstep(i - 1)), 2);
                                rtt{ii} = rt;
                            end
                        end
                        rtmp = cellfun(@(x) x / mean(rtt{1}), rtt, 'uniformoutput', false);
                        tp{i, j} = cellfun(@(x, y) [x; y], tp{i, j}, rtmp, 'uniformoutput', false);
                    end
                    
                    tmp = tp{i, j};
                    subplot(nc - 1, nw, count)
                    mn = cellfun(@(x) nanmean(x, 1), tmp);
                    stdd = cellfun(@(x) nanstd(x, [], 1) / (2 * sqrt(size(x, 1))), tmp);
                    errorbar(mn, stdd)
                    title(ttls{i - 1, j})
                    count = count + 1;
                end
            end
            
        case 'subspace_dimensionality_dims'
            dt = datain{1};
            vals = dt.vals;
            [nc, nw] = size(vals);
            thres = 0.95;
            
            figure(gcf)
            clf
            count = 1;
            for i = 1: nc
                for j = 1: nw
                    tmp = vals{i, j};
                    nm = size(tmp, 1);
                    dimwhwp = NaN(1, nm);
                    dimwhnwp = NaN(1, nm);
                    dimnwhwp = NaN(1, nm);
                    dimnwhnwp = NaN(1, nm);
                    for ii = 1: nm
                        t1 = find(tmp(ii, :, 1) > thres, 1);
                        t2 = find(tmp(ii, :, 2) > thres, 1);
                        t3 = find(tmp(ii, :, 3) > thres, 1);
                        t4 = find(tmp(ii, :, 4) > thres, 1);
                        if ~isempty(t1)
                            dimwhwp(ii) = t1 / 10;
                        end
                        if ~isempty(t2)
                            dimwhnwp(ii) = t2 / 10;
                        end
                        if ~isempty(t3)
                            dimnwhwp(ii) = t3 / 10;
                        end
                        if ~isempty(t4)
                            dimnwhnwp(ii) = t4 / 10;
                        end
                    end
                    
                    subplot(nc, nw, count)
                    mn = [nanmean(dimwhwp), nanmean(dimwhnwp), nanmean(dimnwhwp), nanmean(dimnwhnwp)];
                    st = [nanstd(dimwhwp) / (2 * sqrt(sum(~isnan(dimwhwp)))), nanstd(dimwhnwp) / (2 * sqrt(sum(~isnan(dimwhnwp)))), nanstd(dimnwhwp) / (2 * sqrt(sum(~isnan(dimnwhwp)))), nanstd(dimnwhnwp) / (2 * sqrt(sum(~isnan(dimnwhnwp))))];
                    errorbar(mn, st)
                    ylim([0, 1])
                    count = count + 1;
                end
            end
            
        case 'subspace_dimensionality_ssp_2pair_matrix'
            dt = datain{1};
            ssps = dt.ssps;
            ttls = {'heat full', 'heat no'; 'vf full', 'vf no'};
            [nc, nw, nm] = size(ssps);
            
            figure(gcf)
            clf
            count = 1;
            for i = 1: nc
                for j = 1: nw
                    %%% get tp first: matrix form of all mice %%%
                    %%%% 2-pair %%%%
                    dt2 = zeros(4, 4, nm);
                    for k = 1: nm
                        dt2t = diag(ones(1, 4));
                        for ii = 1: size(ssps{i, j, k}{1}, 1)
                            id1 = str2double(ssps{i, j, k}{1}{ii, 1}{1}(2));
                            id2 = str2double(ssps{i, j, k}{1}{ii, 1}{2}(2));
                            tmp1 = ssps{i, j, k}{1}{ii, 2}{1} / mean(ssps{i, j, k}{1}{ii, 4}{1});
                            dt2t(id1, id2) = tmp1;
                            dt2t(id2, id1) = tmp1;
                        end
                        dt2(:, :, k) = dt2t;
                    end
                    
                    subplot(nc, nw, count)
                    imagesc(mean(dt2, 3))
                    title(ttls{i, j})
                    count = count + 1;
                end
            end

        case 'subspace_dimensionality_ssp_pairs_average'
            dt = datain{1};
            ssps = dt.ssps;
            ttls = {'heat full', 'heat no'; 'vf full', 'vf no'};
            [nc, nw, nm] = size(ssps);
            
            figure(gcf)
            clf
            count = 1;
            for i = 1: nc
                for j = 1: nw
                    %%% get tp first: matrix form of all mice %%%
                    %%%% 2-pair %%%%
                    dt2 = [];
                    for k = 1: nm
                        for ii = 1: size(ssps{i, j, k}{1}, 1)
                            tmp1 = ssps{i, j, k}{1}{ii, 2}{1} / mean(ssps{i, j, k}{1}{ii, 4}{1});
                            dt2 = [dt2, tmp1];
                        end
                    end

                    %%%% 3-pair %%%%
                    dt3 = [];
                    for k = 1: nm
                        for ii = 1: size(ssps{i, j, k}{2}, 1)
                            tmp1 = ssps{i, j, k}{2}{ii, 2}{1} / mean(ssps{i, j, k}{2}{ii, 4}{1});
                            dt3 = [dt3, tmp1];
                        end
                    end

                    %%%% 4-pair %%%%
                    dt4 = [];
                    for k = 1: nm
                        for ii = 1: size(ssps{i, j, k}{3}, 1)
                            tmp1 = ssps{i, j, k}{3}{ii, 2}{1} / mean(ssps{i, j, k}{3}{ii, 4}{1});
                            dt4 = [dt4, tmp1];
                        end
                    end

                    subplot(nc, nw, count)
                    mn = [mean(dt2), mean(dt3), mean(dt4)];
                    st = [std(dt2) / (2 * sqrt(length(dt2))), std(dt3) / (2 * sqrt(length(dt3))), std(dt4) / (2 * sqrt(length(dt4)))];
                    errorbar(mn, st)
                    ylim([0, 1])
                    xlim([0, 4])
                    title(ttls{i, j})
                    count = count + 1;
                end
            end

        case 'demo_manifold_cluster'
            pbf = datain{1};
            mfs = datain{2};
            i = 3;
            j = 1;
            k = 2;
            cl = pbf{i, j, k};
            tp = mfs{i, j, k};
            
            figure(gcf)
            clf
            hold on
            plot3(tp(1, :), tp(2, :), tp(3, :), '.k', 'markersize', 2)
%             for ii = 1: length(cl)
%                 plot3(tp(1, cl{ii}), tp(2, cl{ii}), tp(3, cl{ii}), '.', 'color', clrs(mod(ii - 1, 8) + 1, :), 'markersize', 3)
%             end
            
        case 'toy_manifold_demo_embed'
            f = figure(gcf);
            clf
            [X, labels] = generate_data('swiss', 1000);
            % [X, labels] = generate_data('twinpeaks', 1000);
            ylim([-10, 40])
            dd = squareform(pdist(X));
            thres = prctile(dd(:), 2);
            s = dd < thres;
            s = s - diag(ones(1, 1000));
            s = triu(s);
            [h, w] = find(s);
            idx = [h(:), w(:)];
            idx = idx(randsample(length(h), round(length(h) * 0.5)), :);
            hold on
            for i = 1: size(idx, 1)
                plot3(X(idx(i, :), 1), X(idx(i, :), 2), X(idx(i, :), 3), 'color', [0.7, 0.7, 0.7])
            end
            scatter3(X(:,1), X(:,2), X(:,3), 10, labels, 'filled'); title('Original dataset'), drawnow
            
            f.Renderer = 'Painters';
       
        case 'toy_manifold_demo_flat'
            f = figure(gcf);
            clf
            [X, labels, t] = generate_data('swiss', 1000);
            % [X, labels] = generate_data('twinpeaks', 1000);
            ylim([-10, 40])
            dd = squareform(pdist(X));
            thres = prctile(dd(:), 2);
            s = dd < thres;
            s = s - diag(ones(1, 1000));
            s = triu(s);
            [h, w] = find(s);
            idx = [h(:), w(:)];
            idx = idx(randsample(length(h), round(length(h) * 0.5)), :);
            hold on
            for i = 1: size(idx, 1)
                plot(t(idx(i, :), 1), t(idx(i, :), 2), 'color', [0.7, 0.7, 0.7])
            end
            scatter(t(:,1), t(:,2), 10, labels, 'filled'); title('Original dataset'), drawnow
            
            f.Renderer = 'Painters';

        case 'toy_manifold_demo_learned'
            [X, labels, t] = generate_data('swiss', 1000);
            % [X, labels] = generate_data('twinpeaks', 2000);
            % [X, labels] = generate_data('intersect', 2000);
            scatter3(X(:,1), X(:,2), X(:,3), 5, labels); title('Original dataset'), drawnow
            no_dims = 2;
            % [mappedX, mapping] = compute_mapping(X, 'Isomap', no_dims, 25);
            [mappedX, mapping] = compute_mapping(X, 'Laplacian', no_dims, 10);
            scatter(mappedX(:,1), mappedX(:,2), 5);
            
            [W, L] = compute_geo_laplacian(mappedX, 51);
            [~, rm] = riemann_metric(mappedX, L);
            cover_plot(rm, mappedX, labels)
            
        case 'plot_persistent_barcode'
            pb = datain{1};
            [nc, nw, nm] = size(pb);
            
            f1 = figure(1);
            clf
            f2 = figure(2);
            clf
            count = 1;
            for i = 1: nc
                for j = 1: nw
                    for k = 1: nm
                        dt = pb{i, j, k};
                        tmp = dt.d0;
                        dd = diff(tmp, 1, 2);
                        [~, id] = sort(dd, 'descend');
                        tmp = tmp(id, :);
                        tp = tmp(1: min(50, round(size(tmp, 1) * 0.01)), :);
                        idnan = isnan(tp);
                        tp(idnan) = dt.max;
                        [~, idt] = sort(tp(:, 1), 'descend');
                        tp = tp(idt, :);
                        
                        figure(f1)
                        subplot(6, 6, count)
                        hold on
                        for ii = 1: size(tp, 1)
                            plot(tp(ii, :), 51 - [ii, ii], 'linewidth', 1)
                        end
                        xlim([0, dt.max])
                        
                        tmp = dt.d1;
                        dd = diff(tmp, 1, 2);
                        [~, id] = sort(dd, 'descend');
                        tmp = tmp(id, :);
                        tp = tmp(1: min(50, round(size(tmp, 1) * 0.01)), :);
                        idnan = isnan(tp);
                        tp(idnan) = dt.max;
                        [~, idt] = sort(tp(:, 1), 'descend');
                        tp = tp(idt, :);

                        figure(f2)
                        subplot(6, 6, count)
                        hold on
                        for ii = 1: size(tp, 1)
                            plot(tp(ii, :), 51 - [ii, ii], 'linewidth', 1)
                        end
                        xlim([0, dt.max])
                        count = count + 1;
                    end
                end
            end
    end
end